######################
##### PREPARATION ####
######################
# Clear environment
rm(list=setdiff(ls(), "path_to"))

# Load data
s = full = readRDS(path_to("wide"))
s = s %>% filter(!is.na(code_conseq_simple))


s %>% group_by(conseq, wtp_ratio_class) %>% summarise(
  n = n(),
  wgt = sum(wgt),
  nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
  consumer = sum(consumer_vig),
  dampening = sum(belief_dampening == 1)
) %>%
  mutate(
    freq_main = n / sum(n, na.rm=T),
    freq_wgt = wgt / sum(wgt, na.rm=T),
    freq_nospeeders = nospeeder / sum(nospeeder),
    freq_consumer = consumer / sum(consumer),
    freq_dampening = dampening / sum(dampening),
    wtp_ratio_class=factor(wtp_ratio_class)
  )

######################
#### DERIVE DATA #####
######################
# Define different specifications
specs = list(
  "main" = rep(1, nrow(s)),
  "wgt" = s$wgt,
  "nospeeders" = ifelse(s$Duration__in_seconds_ > quantile(full$Duration__in_seconds_, 0.2), 1, NA),
  "consumer" = ifelse(s$consumer_conseq == 1, 1, NA),
  "dampening" = ifelse(s$belief_dampening == 1, 1, NA)
)

# Anticipated effects
df = list()
for (sp in names(specs)) {
  s$weight = specs[[sp]]
  df[[sp]] = rbind(
    s %>% group_by(conseq) %>% summarise(
      spec = sp,
      freq_conseq = wtd.mean(grepl("conseq", code_conseq_simple) & !grepl("deont", code_conseq_simple), weight),
      freq_both   = wtd.mean(grepl("conseq", code_conseq_simple) & grepl("deont", code_conseq_simple), weight),
      freq_deont  = wtd.mean(!grepl("conseq", code_conseq_simple) & grepl("deont", code_conseq_simple), weight),
      freq_other  = wtd.mean(grepl("other", code_conseq_simple), weight)
    ),
    s %>% summarise(
      spec = sp,
      freq_conseq = wtd.mean(grepl("conseq", code_conseq_simple) & !grepl("deont", code_conseq_simple), weight),
      freq_both   = wtd.mean(grepl("conseq", code_conseq_simple) & grepl("deont", code_conseq_simple), weight),
      freq_deont =  wtd.mean(!grepl("conseq", code_conseq_simple) & grepl("deont", code_conseq_simple), weight),
      freq_other  = wtd.mean(grepl("other", code_conseq_simple), weight)
    ) %>% mutate(conseq = "pooled")
  ) %>% pivot_longer(3:6, values_to = "freq", names_to = "code")
}
df = rbind.fill(df)
df$code = gsub("freq_", "", df$code)
df$code = factor(df$code, levels=c("other", "deont", "both", "conseq"))

# Don't display other fraction
df = df %>% filter(code != "other")

# Order
conseqs = c(
  "pooled" = "All cases",
  "co2" = "CO2 emissions\n(Reduce by 1 ton)",
  "waste" = "Non-recyclable waste\n(Reduce by 100 pounds)",
  "chicken" = "Animal welfare\n(20 certified chicken)",
  "textile" = "Low wages\n(10 fairtrade garments)"
)
df$conseq = factor(df$conseq, levels=names(conseqs))

# Robustness tests
tests = c(
  "main" = 1,
  "wgt" = 2,
  "nospeeders" = 3,
  "consumer" = 4,
  "dampening" = 5
)

######################
#### PLOT ############
######################
plots = list()
plots$fig_explanations_conseq_robust = ggplot(df, aes(x = spec, y = freq, fill = code)) +
  geom_bar(stat = "identity") +
  facet_grid(.~conseq, labeller = as_labeller(conseqs)) +
  xlab("Specification") + ylab("Frequency") +
  coord_cartesian(ylim=c(0, 1)) +
  scale_x_discrete(limits = names(tests), labels=tests) +
  scale_y_continuous(breaks=seq(0, 1, 0.2), labels=paste0(100*seq(0, 1, 0.2), "%")) +
  scale_fill_manual(
    values=c("skyblue", "#B0858D", "red2"),
    labels=c("Action / direct effect  ", "Both  ", "Consequences / total effect"),
    name = "**Explanation**: Respondent cares about ..."
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "rob_con"),
    plot.title = element_text(size=18, face="bold", hjust=0.5),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.title = element_markdown(size=14, hjust=0.5),
    #legend.spacing.y = unit(0, 'cm'),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    strip.text = element_text(size=14, face="bold"),
    axis.title = element_text(size=14, face="bold"),
    axis.text.x = element_markdown(size=14),
    axis.text.y = element_text(size=14),
    legend.text = element_text(size=14)
  ) +
  guides(fill = guide_legend(nrow = 1, reverse = T))

# Save generated graph.
pdf(path_to("figures", "fig_explanations_conseq_robust.pdf"), width=11, height=5.5)
print(plots$fig_explanations_conseq_robust)
dev.off()

